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ABSTRACT 

We exploit a generic instability in the integration of steady, sonic, near-isothermal 
flows to find the complete transition diagram for recombination fronts (for a model 
system of equations). The instability requires the integration of the flow equations 
for speeds between the isothermal and adiabatic sound speeds to be performed with 
particular care. As a result of this, the previous work of Newman & Axford on the 
structure of recombination fronts neglected an important class of solution, that of 
transonic fronts; our method is readily extensible to a more complete treatment of the 
ionization structure. Future papers will apply these results in models of the structure 
of ultracompact Hn regions. 
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1 INTRODUCTION 

The hydrodynamics of gas pho toionized by a central star 
has been widely s tudied (e.g. Kahn 1954; Axford 1961 



Goldsworthy 1961). The gas is almost completely ionized 



close to the central source. The incident radiation field de- 
creases with radius as a result of geometrical divergence and 
because atomic recombination always results in a small pop- 
ulation of neutral atoms, which have finite opacity to the 
ionizing radiation. At some radius the gas becomes opti- 
cally thick to the ionizing radiation, and beyond this radius 
it is predominantly neutral. 

The transition between ionized and neutral gas can be 
treated as a quasi-steady, plane-parallel front when it takes 
place over a small distance, compared with the radius of 
the region, and can be i ncluded in models of the global flow 
as a discontinuity (e.g. psterbrock 1989] ). This is, in gen- 
eral, the appropriate case for ionization fronts (IF), since 
the opacity of the gas to ionizing radiation increases rapidly 
as a function of decreasing ionization. In this paper, we will 
distinguish recombination fronts (RF) from IF as the case 
when the net flow of gas is from high ionization towards low 
ionization. RF can be broadened significantly by the advec- 
tion of ionized gas, whereas IF are generally only as thick as 
the (small) photon mean free path in neutral gas. 

Newman & Axford (1968 - henceforth NA) studied the 
possible RF transitions in the context of stellar winds. They 
found that thin fronts were unlikely to occur in the winds of 
solar-type stars but might be important in the ejecta of plan- 
etary nebulae . In recent studies of ultracompact H n regions 
(bvson 1994 jWilliams & Dyson 1994 |Dyson, Williams & 
Redman 1995 ), we have discussed the likely importance of 
(thin) RF in the flow structures. In particular, the likely co- 
incidence of RF and sonic transitions raised some questions 



with regard to the work of NA, which we address in this 
paper. 

The significance of the small class of steady solutions 
which pass through sonic transitions at the RF is perhaps 
not immediately apparent when only conditions close to the 
front are considered. However, such solutions become impor- 
tant when the front is included in the global solution of an 
initially subsonic divergent flow. Here, the ballistic character 
of the supersonic solutions allows such 'winds' to satisfy a 
wide range of external pressure boundary conditions, by the 
inclusion of a shock at some radius, for external pressures 
far lower than required if the solutions are subsonic at all 



radii ( Parker 1958| ). 

As part of this study, we also address the apparent para- 
dox that for conditions arbitrarily close to isothermal, the 
flow equations are singular at the adiabatic sound speed, 
whereas in the isothermal limit they are singular at the 
isothermal sound speed. 

In the following sections, we give simple model equa- 
tions for steady ionized gas flow, and derive a second order 
system which describes the structure of plane IF and RF. 
A reorganisation of these equations suggests that for a nar- 
row range of flow speeds in the region where the flow is 
nearly isothermal, a formal instability will prevent direct in- 
tegration of the Mach number equation in the downstream 
direction; this is illustrated by a yet simpler model system. 
With this in mind, we return to RF, and find a range of tran- 
sonic solutions which were inadvertently missed by NA. In 
conclusion, we briefly discuss the relevance of these results 
to models of ultra compact Hn regions, and note possible 
implications of our study for the details of the evolution of 
conventional IF in H n regions at the time of shock genera- 
tion. 
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2 BASIC EQUATIONS 

Differential equations for the structure of ionization and re- 



combination fronts are given by NA (cf. also Frank, Noriega- 
Crespo & Balick 1992). Terms which account for photoion- 
ization and recombination, the radiative transfer of ioniz- 
ing radiation and the heating and cooling of the gas are 
added to the basic equations of hydrodynamics. The pro- 
cesses included are ionization of atomic hydrogen, proton 
and electron recombination in the on-the-spot approxima- 
tion, an d an approximation to forbidden line cooling (Lei). 
Mason ( 1977 ) gives far more detailed models of the internal 
structure of several types of IF, showing, for example, the 
ionization structure of O, N and Ne. For the present paper, 
however, the approximations of NA will suffice. 

The following equations give a reasonable approxima- 
tion to the structure of RF , in the frame in which the front 



is steady, (cf. [Axford 1961 ; NA ) 
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where p is the mass-density in the gas, v is the velocity, p is 
the pressure and x is the fractional ionization. The equation 
of state, 



p = (1 + x)nkT, 



(2) 



relates these quantities to n — p/mu, the nucleon number 
density, and the gas temperature T. The photoionization 
cross-section of hydrogen at the Lyman limit is taken to 
be a = 6.3xl0 _18 cm 2 per atom, and /3, the case B re- 
combination coefficient, is /?oT _3//4 , where we use the value 
/3o = 3x lCT 10 cm 3 s- 1 K 3/4 from NA. 

The internal energy per unit volume is e = p/{"/ — 1); 
the enthalpy per unit volume is w — e + p. J = Jf is the 
outward flux of ionizing photons. 

Q is an effective heating rate, given by 



pQ = fcT*n(aJ)(l - x) 



7-1 



PkTn 2 x 2 — n 2 x 2 L c i, (3) 



where T* is the effective stellar temperature (which we take 
to be 40 000 K). The final term in equat ion (|^) is a fit to 
forbidden line cooling ( Goldsworthy 1961), where 



L c i — 



L ci0 (T- 4000 K) 2 




for T > 4000 K 
otherwise 



(4) 



and L ci0 = 9.7xl0~ 32 erg cm 3 s _1 K -2 . This fit is poor for 
gas temperatures above 10 K, but since temperatures as 
high as this are not found in the fronts that we are studying, 
and since we wish to compare our results to those of NA, we 
defer a more accurate treatment of forbidden line cooling to 
the future. 

Boundary conditions (at small r) are of high ionization 
(a; — > 1), high ionizing flux (J — > oo), equilibrium tempera- 
ture (Q — * 0) and an arbitrary initial velocity (v > 0). These 
are, however, not sufficient to fully determine the flow in all 
cases, as will be described below. 



3 STEADY STATE SOLUTIONS 

For a thin front, the mass flux and momentum equations ( [la] , 
|rpj ) can be integrated immediately. The remaining ordinary 
differential equations which describe the front are 
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where c is the adiabatic sound speed, c 2 = jp/p, m is the 
adiabatic Mach number, and j = J/nv is the dimensionless 
ionizing flux. The integrals of the mass and momentum flux 
equations, $ and II respectively, are 
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Hence n/<I> = (1 + jm 2 )c/jm. In this system of equations, 
the radial distance scales simply with density, and so we can 
study two coupled ODEs for dm/dj and dx/dj which are 
equivalent to the equations studied by NA. 



3.1 A numerical instability of near-isothermal 
sonic flows 

Using equations ^ to relate m to c, we can rearrange equa- 
tion (pal) into the suggestive form: 



dln(c 2 
dr 



(7) 



(cf. also equation (15) of Williams, Hartquist & Dyson, 
1995). If the level of ionization is weakly varying, then 
c 2 oc T. For substantially subsonic flows (m 1) and for 
supersonic flows (m > 1) the solution will tend to relax 
to a stable equilibrium as we integrate in the downstream 
direction. For our system of equations, the flow reaches equi- 
librium at a temperature T ~ 7800 K (cf. NA). 

However, when the flow velocity is in the range be- 
tween the isothermal sound speed (i.e. m = lA/7) and 
the adiabatic sound speed (m = 1), the first factor on 
the r.h.s. of equation (W) is negative: thermally stable equi- 
libria appear to become unstable (when the flow is inte- 
grated downstream) and vice versa (so long as the dynamical 
time, (du/dr)" 1 , is much longer than the thermal timescale, 
~c 2 /Q). 

This behaviour is readily understood in terms of the for- 
mation of weak 'isothermal shocks' in the flow. If the flow 
were strictly isothermal, discontinuous shocks could form 
with upstream Mach numbers in the unstable range (cf. Fig- 
ure ^J) . For finite cooling, discontinuous solutions no longer 
exist. The instability of the equilibrium solution neverthe- 
less allows the flow to decrease its Mach number through 
a resolved transition at any point, and since the flow prop- 
erties upstream and downstream of this transition satisfy 
the isothermal Rankine-Hugoniot relations, we refer to these 
transitions as resolved isothermal shocks. 

We illustrate this in Figure [l] for the simpler case of a 
steady plane flow with no ionization, using a simple heat- 
ing function Q oc (T c — T) in equation (pa|) , where T c is the 
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(a) (b) 




Figure 1. Steady solutions of flows with heating Q oc T c — T for (a) moo = 1.2 and (b) TTioc — 0.9. The flow direction is towards 
increasing r. The bold, horizontal lines are the equilibrium solutions. Dotted lines mark the isothermal and adiabatic sound speeds. 



equilibrium temperature. The plot shows rest frame solu- 
tions, with the integration performed downstream, towards 
increasing r. The solutions are independent of r, but a range 
of solutions are shown arbitrarily offset in this direction to 
indicate the morphology of solutions expected in more com- 
plicated situations. 

Note that in this model, Q is independent of p. This 
cooling law is stable to equilibrium perturbations, and also 
marginally stable to the instab i lities of algebraic power 



istence of rapidly divergent solutions for flows between the 
isothermal and adiabatic sound speeds, we c onclude that the 
solutions found by Newman & Axford ( NA ) in this range of 



within cooling zones (Field 1965; Fall & Rees 1985). Since 
the peak temperatures found within the isothermal shocks 
are less than 3 per cent above the equilibrium temperature, 
we can safely ignore this marginal instability; in any case, 
the results are similar for Q oc (T™ — T n ) with n > 1. 

We first show a case for which both the sub- and the 
supersonic solutions to the flow equations integrate stably, 
with moo = 1.2 (or, equivalently, m,oo = 0.5). In Figure |l|a, 
the bold, horizontal lines are the equilibrium solutions; thin- 
ner solid lines show solutions integrated downstream (to the 
right) from an initial non-equilibrium state. The isothermal 
and adiabatic sound speeds are shown as dotted lines. As we 
integrate downstream, sub- (super-) sonic initial conditions 
relax to the sub- (super-) sonic equilibrium solution. 

Figure |l|b is for moo = 0.9 (or, equivalently, moo — 
0.67). Here, we see that the equilibrium solution m = 0.9 
is apparently unstable as we integrate downstream (left to 
right), with solutions diverging to large m and towards the 
m ~ 0.67 equilibrium. These latter solutions are, in effect, 
weak shocks. If we consider a Galilean transformation of the 
system, it is apparent that the instability found is entirely 
a result of the steady-state assumption: m = 0.9 is a per- 
fectly reasonable equilibrium solution, as is confirmed by 
time-dependent simulations of such flows. For some down- 
stream conditions, however, it will be necessary for a weak 
shock to form in the flow at some (in the present case, ar- 
bitrary) position. The instability demonstrated here allows 
this to occur in continuous solutions, just as discontinuous 
adiabatic Rankine-Hugoniot shocks can form at arbitrary 
positions in flows initially above the adiabatic sound speed. 

Since downstream integration of these ODEs for ini- 
tial conditions perturbed from equilibrium confirms the ex- 



upstream Mach numbers inadvertently included these weak 
isothermal shocks. In the next section, we present the addi- 
tional transonic solutions this led them to disregard. 



3.2 Structure of recombination fronts 

The presence of the instability discussed in the previous 
section in the downstream integration of near isothermal 
flows implies that, when the upstream flow is between the 
isothermal and adiabatic sound speeds, the initial conditions 
m — > moo, x — > 1 at j — > oo used by Newman & Axford do 
not fully specify the flow solutions. In this section, we use 
this freedom to find an important class of RF overlooked by 
NA - transonic fronts. In what follows, we return to the full 
model equations (|^a-c) for the structure of RF. 

In looking for transonic solutions, where Q = and 
v(dQ/dr) < at m = 1 (so the flow can pass smoothly 
through the adiabatic sound speed), we have formulated a 
two-point boundary problem (with conditions applied both 
far upstream in the Hn flow and at the sonic point). We 
have indicated that for initial conditions l/y^y < m < 1, 
the Mach equation is unstable when integrated downstream. 
However, the ionization equation is highly unstable when 
integrated in the opposite, upstream direction. 

The vigorousness of these instabilities precludes the use 
of shooting methods. We therefore employed an iterative 
procedure to find the flow solution between j — > oo and 
the adiabatic sonic point (outside this region, downstream 
integration proceeds smoothly). An approximation to the 
ionization structure was determined by downstream inte- 
gration. Using this approximation to the ionization struc- 
ture, the equation for the Mach number was integrated up- 
stream. This gave an updated model for the velocity and 
sound speed distributions, which were then used to derive 
an updated ionization structure: this procedure was iterated 
to convergence. 

The integratio ns were performed with a stiff ODE in- 
tegration package ( Cohen fc Hindmarsh 1994 ), using cubic 
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Figure 3. Example solutions for moo = 0.85. The enthalpy (h, in units of c^), sound speed (c in units of Coo) and adiabatic Mach 
number (m) are plotted against the flow velocity [u in units of Coo). Vertical dotted lines mark the maximum of c, where the flow is at 
the isothermal sound speed, and the maximum of h, where the flow is at the adiabatic sound speed; horizontal dotted lines mark the 
values of m at these speeds. In case (a) an 'isothermal shock' is formed ahead of the sonic point (the flow travels A — > B then passes 
through the isothermal shock B — > A — > B' and then recombines and cools to reach C in the end). Case (b) includes a sonic point (S) — 
the solution with no shock is A — > S — » B — > C, but there are also solutions such as A— >S^B^D^E which include an adiabatic 
shock (B — > D, shown dashed). 



spline interpolation of the results. The singular equations 
for dm/dj and dx/dj were split into non-singular equations 
in a dummy independent variable s (i.e. for dm/ds, dx/ds 
and dj/ds). 

In Figure ^| we show one such transonic solution, for 
upstream Mach number rrtoo = 0.85. The solid bold curve 
is the Mach number of the transonic solution, plotted as 
a function of the ionization fraction, x. Also shown in this 
figure are the distributions of dimensionless ionizing flux j, 
sound speed c, temperature T, velocity v and stagnation 
enthalpy, h — c 2 /(-f — 1) + f 2 /2, through the front. Note 
that all these curves are smooth, and that h has a maximum 
value when m = 1, as expect (cf. Figure |^). 

In Figure ^, we show the various classes of solution per- 
turbed from the transonic case. The three curves on the plots 
show the adiabatic Mach number, m, sound speed, c, and 



enthalpy ft as a function of flow velocity, v. In Figure ga, a 
solution with a weak isothermal shock is shown. From the 
initial point A, the flow cools and accelerates. However, an 
isothermal shock is triggered at point B, and the flow then 
heats again and slows. As the shock is resolved, the solution 
remains on the integral curves at all times - the part of the 
curve between points B and B' might, however, be consid- 
ered to represent the isothermal shock. Beyond the shock, 
in the region B' — > C, the flow decelerates and cools further 
as the gas recombines. 

Figure |^b illustrates a transonic solution, A — > S — > 
B — > C. Here, the flow has been chosen to accelerate 
smoothly through sonic point S, at the maximum of the en- 
thalpy h. Beyond the sonic point, the flow accelerates, cools 
and recombines to reach point C when it has become fully 
neutral. Also illustrated is the solution which results when 
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Figure 4. Steady solutions for an RF with Mach number m^c = 0.85 far upstream (at high ionizing flux, to the left of this plot) in 
the ionized gas. The bold curve is the transonic solution. A sample of the solutions divergent from the transonic case are shown for 
initial points upstream of the sonic point. Downstream of the sonic point, solutions are shown, perturbed by an initial steady adiabatic 
shock (shown dashed). Horizontal dotted lines mark the isothermal and adiabatic sound speeds. For the subsonic downstream solutions, 
final temperatures vary between 6500 K (for flows which shock far upstream in the H II region) and 14000 K (where the gas shocks far 
downstream in the Hi region). 



an adiabatic shock occurs at point B, taking the flow dis- 
continuously to D. Beyond D, the flow behaves in a similar 
fashion to that illustrated for the isothermal shock case in 
Figure |§|a, except that it does not reach quite such a low 
temperature as it approaches infinity. 

In Figure ^ we plot the flow Mach number, m, for a 
wide range of these solutions against the dimensionless ion- 
izing flux, j. Upstream of the sonic point, both the phys- 
ical solution (which continues subsonically to j — > 0) and 
the unphysical solution (which tends to m = 1 with infi- 
nite gradient) are shown. At large j, these solutions diverge 
very sharply, underlining the difficulty in finding transonic 
solutions by shooting methods. Beyond the sonic point, so- 
lutions are shown for initial conditions found downstream 
of steady adiabatic shocks (shown dashed) which are now 
allowed. These shocked solutions form a one-parameter fam- 
ily in combination with the physical solutions which diverge 
from the transonic solution upstream of the sonic point: as 
the stepping-off point of these solutions moves progressively 
downstream, the maximum in Mach number of the diver- 
gent solutions becomes a cusp at the sonic point, and then 
a discontinuity (i.e. shock). A range of final Mach numbers 
results from this range of solutions. 



In Figure |5| we show the transition diagram from Mach 
number 777,2 in the Hn flow (at the equilibrium tempera- 
ture) far upstream to Mach number mi in the downstream 
Hi flow. In this diagram, we plot both the results we ob- 
tained following Newman & Axford, and the extra solutions 
found in this work. The short segment of the upper right 
hand branch of solutions shown bold, between dotted hor- 
izontal lines at the isothermal and adiabatic sound speeds, 
give the transonic solutions. Additional regions of subsonic 
solutions to the top left result both from shocks downstream 
of the sonic point in transonic solutions, and from the use of 
isothermal, rather than adiabatic, shock relations for shocks 
forming in the H n region. 

To interpret this diagram, it is useful to remember that 
because the upstream flow is assumed to be at the fixed 
equilibrium temperature, the curves plotted are for a family 
of ODEs, parameterised by II/<E>. Isothermal shock condi- 
tions m+m_ = 1/7 relate the values of m,2 where the ODEs 
are identical. The sudden appearance of extra solutions with 
mi > 1 at 7712 = 1/ ^7 i s allowed because this is the singular 
value where 11/ $ is a minimum. 
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Figure 2. Transonic solution for an RF with Mach number 
rrica = 0.85 far upstream in the ionized gas. The Mach num- 
ber m, dimensionless flux j, sound speed c, temperature T, flow 
velocity v and specific enthalpy h are plotted against ionization 
fraction x. This solution here has j Bon = 0.51 and x son = 0.936, 
and the Mach number of the H I flow far downstream is mi = 2.56, 
corresponding to a temperature there of 4800 K. 




Figure 5. Transition diagram for RF (cf. Figure 5 of NA). Here 
ni2 is the adiabatic Mach number far upstream in the H II flow, 
mi is the adiabatic Mach number far downstream in the H I flow. 
The dashed line is mi = ni2\ dotted lines mark where the up- 
stream flow is at the isothermal sound speed, m,2 = X/^y, and 
at the adiabatic sound speed m2 = 1- The narrow lines and dark 
gray area are the same as those in Newman & Axford (although 
note that they use isothermal, rather than adiabatic, Mach num- 
bers). The bold lines and lighter shaded area are also allowed. 
In particular, note the short section of the mi > 1 curve with 
between the horizontal dotted lines: this section is the locus of 
transonic RF. 



4 CONCLUSIONS 

In this paper, we have comp leted the solution space found 
by Newman & Axford (NA) for recombination front flows. 
Careful treatment of the numerical integration for flow 
speeds between the isothermal and adiabatic sound speeds 
reveals an extra class of transonic solutions overlooked by 
these authors. 

Close to the sonic point, the form of these solutions is 
similar to those found for any steady transonic flow, for in- 



tions is at the adiabatic sound speed, whereas in a strictly 
isothermal flow it is at the isothermal sound speed. We find 
that, where thermal timescales are short in an accelerating 
flow, a 'wide' range of solutions near the adiabatic sound 
speed derive from a 'few' at the isothermal sound speed: it 
is apparent that in the limit, these 'few' solutions can be- 
come the singular isothermal transonic solution. 

In the theory of multifluid MHD shocks, smo oth tran- 
sonic solutions are also found for strong cooling (Chernofl 



can in j;eneral be found by integrating away from the sonic 



stance i Parker (|l958|) wind. The steady, transonic solution 1987 ; Flower & Pineau des Forets 1987; Robcrgc & Drains 



199C). In these shocks, the solution adapts to pass through 



point. However, inclusion of ionization equilibrium means 
that we must treat a two-point boundary value problem. 
The strength of the instability discussed above adds signif- 
icantly to the usual difficulties integrating the Mach num- 
ber equation towards a nodal solution, rendering shooting 
methods impractical and allowing transonic solutions to be 
easily overlooked. We overcome this difficulty with an itera- 
tive scheme, integrating separate equations in the upstream 
and downstream directions. 

The instability in downstream integration adds a free 
parameter to the solution, resolving the conceptual difficulty 
of how, as the initial Mach number of our solution increases 
through isothermal Mach 1, the flow solution changes from 
an initial- value to a two-point boundary problem. Moreover, 
these solutions give some insight into the apparent paradox 
that, for finite cooling, the singularity in steady flow equa- 



a downstream (saddle type) sonic point using the freedom 
of manoeuvre allowed by an upstream viscous subshock (J- 
type) or nodal sonic point (C*-type), rather than as dis- 
cussed here. 

There has been little discussion of RF since the work 
of NA, since in the winds of main sequence stars they will 
occur, if at all, at large radii, and thus at densities and 
temperatures so low that they will be broadened beyond 
recognition by advection. The interstellar medium in the 
vicinity of strongly ionizing stars does not generally reach 
an equilibrium state in the stellar lifetime - an IF is still 
being driven out into the surrounding medium when the 
star dies. However, in recent work we have discussed how a 
near-equilibrium state may found in the youn gest of these 
stars, observed as ultra compact Hii regions ( Dv^on^lJWj ; 
Williams fc Dyson 199^ ; |Dyson, Williams fc Redman 1995| ). 
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Mass loading in the near vicinity of the young stars traps the Williams R.J.R., Dyson J.E., 1994, MNRAS, 270, L52 

IF at a scale of 0.1 pc, turning it into a RF through which 

gas passes at high density. In some objects, the mass-loaded This paper has been produced using the Royal Astronomical 
stellar wind may remain supersonic at all radii - however, Society/Blackwell Science fflrjX style file, 
if the mass loading is strong enough the flow will shock. 
Detailed hydrodynamical simulations (Williams, Redman & 
Dyson, in preparation) confirm that transonic RF are fre- 
quently a component of the equilibrium flow. Such regions 
may be a significant source for partially ionized gas in the 
interstellar medium. 

We note briefly that the discussion here also applies to 
the structure of H n regi on IF as they undergo the transition 
from R-type to D-type ( Mason 1980). A shock first forms, as 
a (quasi-) steady part of the IF structure, at the minimum 
of the adiabatic Mach number in the R-type front when this 
reaches unity; initially, the post-shock flow is in the unstable 
range. In this region, however, the heating and cooling in the 
transonic flow are rather weak. More significantly, as the re- 
sultant strong D-type front weakens, the velocity relative to 
the front downstream in the H n region will eventually reach 
the unstable range. In this case, there is a significant chance 
that solutions will be missed in steady-state integrations. 
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